clear all
set more off
cd /Users/zimaoxiao/Dropbox/CC_Yield_Predict/replication_package/  // Navigate to replication folder on your own machine
global weather "timeInt0_3 timeInt3_6 timeInt6_9 timeInt9_12 timeInt12_15 timeInt15_18 timeInt18_21 timeInt21_24 timeInt24_27 timeInt27_30 timeInt30_33 timeInt33_36 timeInt36_39 prec prec2"

qui {
* ============================
* 	  Linear trends spec.
* ============================
use if inrange(year,1950,2000) using data/regression/reg_2020_prediction_3Cstep, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
	
* regression
reghdfe lncornyield $weather state_fips#c.trend [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_lin
parmest, label saving(`"`corn_st_lin'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)

* extract temperature coefficients		
use parm estimate using `corn_st_lin', clear
keep if _n<=13
forvalues i = 0(3)36 {
	local j = `i' + 3
	gen timeInt`i'_`j' = estimate if parm == "timeInt`i'_`j'"
	sum timeInt`i'_`j'
	local k = r(mean)
	replace timeInt`i'_`j' = `k'
}
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace

* extract county linear trends coefficients		
use parm estimate using `corn_st_lin', clear		
drop if _n<16
gen quad = ustrright(parm,1)
drop if quad != "d"
drop quad
destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
drop parm			
rename estimate stlin
tempfile cornlintrend
save `cornlintrend', replace
						
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2020_prediction_3Cstep, clear
bys county_fips: keep if _n==1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
cross using `corntempcoef'
	
* calculate climate change impacts by county
bys county_fips: gen ccimpact_corn = (CC_timeInt0_3*timeInt0_3) + (CC_timeInt3_6*timeInt3_6) + (CC_timeInt6_9*timeInt6_9) + (CC_timeInt9_12*timeInt9_12) + (CC_timeInt12_15*timeInt12_15) + (CC_timeInt15_18*timeInt15_18) + (CC_timeInt18_21*timeInt18_21) + (CC_timeInt21_24*timeInt21_24) + (CC_timeInt24_27*timeInt24_27) + (CC_timeInt27_30*timeInt27_30) + (CC_timeInt30_33*timeInt30_33) + (CC_timeInt33_36*timeInt33_36) + (CC_timeInt36_39*timeInt36_39)
				
* calculate overall impacts by county (climate change + tech effects)
* 1950 as the 1st year, t=1; 2000 as the 51st year, t=51; 2020 as the 71st year, t=71
bys county_fips: gen trends_corn = (CC_timeInt0_3*timeInt0_3) + (CC_timeInt3_6*timeInt3_6) + (CC_timeInt6_9*timeInt6_9) + (CC_timeInt9_12*timeInt9_12) + (CC_timeInt12_15*timeInt12_15) + (CC_timeInt15_18*timeInt15_18) + (CC_timeInt18_21*timeInt18_21) + (CC_timeInt21_24*timeInt21_24) + (CC_timeInt24_27*timeInt24_27) + (CC_timeInt27_30*timeInt27_30) + (CC_timeInt30_33*timeInt30_33) + (CC_timeInt33_36*timeInt33_36) + (CC_timeInt36_39*timeInt36_39) + (stlin*20)	
keep county_fips ccimpact_corn trends_corn 
tempfile corn_cccoef
save `corn_cccoef', replace
	
* merge in regression coefficients, calculate predicted 2020 yields
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef', keep(3) noreport nogenerate
	
* predicted 2020 yields: climate change
bys county_fips: gen cc2020_lin = cornyield_2000*(1+ccimpact_corn)

* predicted 2020 yields: climate change + tech effects
bys county_fips: gen trends2020_lin = cornyield_2000*(1+trends_corn)

* Table S2: rows (1), (2), (4) and (6)		
rename (cornyield_2000 cornyield_2020) (yield2000 yield2020)
keep county_fips yield2020 yield2000 cc2020_lin trends2020_lin
tempfile baseline_2020_lin
save `baseline_2020_lin', replace

* ============================
* 	 Quadratic trends spec.
* ============================
use if inrange(year,1950,2000) using data/regression/reg_2020_prediction_3Cstep, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
				
* regression
reghdfe lncornyield $weather state_fips#c.trend state_fips#c.trend2 [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_quad
parmest, label saving(`"`corn_st_quad'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)
											
* extract temperature coefficients		
use parm estimate using `corn_st_quad', clear
keep if _n<=13
forvalues i = 0(3)36 {
	local j = `i' + 3
	gen timeInt`i'_`j' = estimate if parm == "timeInt`i'_`j'"
	sum timeInt`i'_`j'
	local k = r(mean)
	replace timeInt`i'_`j' = `k'
}
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace
						
* extract county linear trends coefficients		
use parm estimate using `corn_st_quad', clear
drop if _n<16
gen quad = ustrright(parm,1)
preserve
	* first-order
	drop if quad != "d"
	drop quad
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	drop parm			
	rename estimate stlin
	tempfile cornlintrend
	save `cornlintrend', replace
restore
	* second-order
	keep if quad == "2"
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	gen state_fips_1 = floor(state_fips/10)
	replace state_fips = state_fips_1
	drop state_fips_1
	drop parm quad
	rename estimate stquad
	tempfile cornquadtrend
	save `cornquadtrend', replace
	
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2020_prediction_3Cstep, clear
bys county_fips: keep if _n==1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
merge m:1 state_fips using `cornquadtrend', assert(3) keep(3) noreport nogenerate
cross using `corntempcoef'
	
* calculate climate change impacts by county
bys county_fips: gen ccimpact_corn = (CC_timeInt0_3*timeInt0_3) + (CC_timeInt3_6*timeInt3_6) + (CC_timeInt6_9*timeInt6_9) + (CC_timeInt9_12*timeInt9_12) + (CC_timeInt12_15*timeInt12_15) + (CC_timeInt15_18*timeInt15_18) + (CC_timeInt18_21*timeInt18_21) + (CC_timeInt21_24*timeInt21_24) + (CC_timeInt24_27*timeInt24_27) + (CC_timeInt27_30*timeInt27_30) + (CC_timeInt30_33*timeInt30_33) + (CC_timeInt33_36*timeInt33_36) + (CC_timeInt36_39*timeInt36_39)
				
* calculate overall impacts by county (climate change + tech effects)
* 1950 as the 1st year, t=1; 2000 as the 51st year, t=51; 2020 as the 71st year, t=71
bys county_fips: gen trends_corn = (CC_timeInt0_3*timeInt0_3) + (CC_timeInt3_6*timeInt3_6) + (CC_timeInt6_9*timeInt6_9) + (CC_timeInt9_12*timeInt9_12) + (CC_timeInt12_15*timeInt12_15) + (CC_timeInt15_18*timeInt15_18) + (CC_timeInt18_21*timeInt18_21) + (CC_timeInt21_24*timeInt21_24) + (CC_timeInt24_27*timeInt24_27) + (CC_timeInt27_30*timeInt27_30) + (CC_timeInt30_33*timeInt30_33) + (CC_timeInt33_36*timeInt33_36) + (CC_timeInt36_39*timeInt36_39) + (stlin*20) + (stquad*2440)			
keep county_fips ccimpact_corn trends_corn 
tempfile corn_cccoef
save `corn_cccoef', replace
	
* merge in regression coefficients, calculate predicted 2020 yields
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef', keep(3) noreport nogenerate
				
* predicted 2020 yields: climate change
bys county_fips: gen cc2020_quad = cornyield_2000*(1+ccimpact_corn)
				
* predicted 2020 yields: climate change + tech effects
bys county_fips: gen trends2020_quad = cornyield_2000*(1+trends_corn)

* Table S2: rows (3) and (5)		
rename (cornyield_2000 cornyield_2020) (yield2000 yield2020)
keep county_fips cc2020_quad trends2020_quad
tempfile baseline_2020_quad
save `baseline_2020_quad', replace

* ============================
* 	         Output
* ============================
use `baseline_2020_lin', clear
merge 1:1 county_fips using `baseline_2020_quad', assert(3) keep(3) nogenerate
noi logout, save(output/TableS2) word dec(1) replace: tabstat yield2000 cc2020_lin cc2020_quad trends2020_lin trends2020_quad yield2020, stat(p10 p25 p50 p75 p90) format(%10.4f) c(s)
save data/figure/baseline_2020, replace // for Figure 2

}

*** EOF
